#### installing packages if not installed ####

list.of.packages <- c('gridExtra',
                      'dplyr',
                      'ggthemes',
                      'rdrobust',
                      'data.table',
                      'estimatr',
                      'tidyverse',
                      'haven',
                      'readxl')
new.packages <- 
  list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

#### attach packages ####

suppressPackageStartupMessages(
  
  {
    
    library(gridExtra)
    library(dplyr)
    library(ggthemes)
    library(rdrobust)
    library(data.table)
    library(estimatr)
    library(tidyverse)
    library(haven)
    library(readxl)
    
  }
  
)


#### cleaning data ####

load(file = 'data/nyc_cleaned.RData')

df_stops_nyc_ts = df_stops_nyc %>%
  mutate(por = ifelse(date >= as.Date('2020-06-15'), 1, 0)) %>% 
  mutate(porrv = date - as.Date("2020-06-15")) %>% 
  mutate(blm = ifelse(date >= as.Date('2020-05-28'), 1, 0)) %>% 
  mutate(blmrv = date - as.Date("2020-05-28")) %>% 
  group_by(date) %>% 
  summarize(count = sum(count),
            blm = mean(blm),
            blmrv = mean(blmrv),
            por = mean(por),
            porrv = mean(porrv),
            arrest = mean(arrest),
            hit = mean(hit),
            r_blk = sum(r_blk),
            r_wht = sum(r_wht),
            r_lat = sum(r_lat)) %>% 
  mutate(rate_blk = (r_blk / (.228 * 8175133)) * 10000,
         rate_wht = (r_wht / (.333 * 8175133)) * 10000,
         rate_lat = (r_lat / (.286 * 8175133)) * 10000) %>% 
  mutate(rr_bw = rate_blk / (rate_wht+mean(rate_wht, na.rm = TRUE)),
         rr_lw = rate_lat / (rate_wht+mean(rate_wht, na.rm = TRUE)))

#### descriptive justification #### 

descjus1 = df_stops_nyc %>% 
  group_by(date) %>% 
  summarize(count = sum(count),
            blm = mean(blm),
            blmrv = mean(blmrv)) %>% 
  filter(date >= as.Date("2020-03-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = count, group = blm),
              col = 'black',
              se = FALSE,
              size = .6) + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2) + 
  annotate('text',
           x = as.Date("2020-07-01"),
           y = 40,
           label = "BLM Protest\nOnset\n(2020-05-28)",
           family = 'serif',
           size = 3) + 
  labs(x = "Date", y = "Stops",
       title = "A. BLM Protest Onset") + 
  theme_tufte()

mod2ann = 
  with(df_stops_nyc_ts, rdrobust(x = porrv, y = count, p = 1, kernel = 'uni'))

descjus2 = df_stops_nyc %>% 
  group_by(date) %>% 
  mutate(blm = ifelse(date >= as.Date('2020-06-15'), 1, 0)) %>% 
  summarize(count = sum(count),
            blm = mean(blm),
            blmrv = mean(blmrv)) %>% 
  filter(date >= as.Date("2020-03-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_point(aes(x = date, y = count),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = count, group = blm),
              col = 'black',
              se = FALSE,
              size = .6) + 
  geom_vline(xintercept = as.Date("2020-06-15"),
             linetype = 2) + 
  annotate('text',
           x = as.Date("2020-07-07"),
           y = 50,
           label =   paste0('RDiT Est: ', round(mod2ann$coef[3], 2), "\nSE: ",
                            round(mod2ann$se[3], 2),
                            '\np < 0.001'),
           family = 'serif',
           size = 2.5) + 
  annotate('text',
           x = as.Date("2020-07-07"),
           y = 33,
           label = "Plainclothes Officer\nDissolution\n(2020-06-15)",
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Stops",
       title = "B. Plainclothes Officer Dissolution") + 
  theme_tufte()

  
descjusgrob = arrangeGrob(descjus1, descjus2, ncol = 2)
ggsave(plot = descjusgrob, filename = 'pics/descjus.jpeg',
       dpi = 600, width = 8, height = 2.5)

#### main results: policing quality ####

mainres_stops = 
  with(df_stops_nyc_ts, rdrobust(x = porrv, y = count, p = 1, kernel = 'uni'))
mainres_arr = 
  with(df_stops_nyc_ts, rdrobust(x = porrv, y = arrest, p = 1, kernel = 'uni'))
mainres_hit = 
  with(df_stops_nyc_ts, rdrobust(x = porrv, y = hit, p = 1, kernel = 'uni'))
mainres_rrbw = 
  with(df_stops_nyc_ts, rdrobust(x = porrv, y = rr_bw, p = 1, kernel = 'uni'))
mainres_rrlw = 
  with(df_stops_nyc_ts, rdrobust(x = porrv, y = rr_lw, p = 1, kernel = 'uni'))

resny1 = 
  df_stops_nyc_ts %>% 
  filter(date >= as.Date("2020-04-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-06-15"),
             linetype = 2) + 
  geom_point(aes(x = date, y = hit),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = hit, group = por),
              col = 'black',
              se = FALSE,
              size = .6) + 
  annotate('text',
           x = as.Date("2020-07-07"),
           y = .27,
           label = paste0('RDiT Est: ',
                          round(mainres_hit$coef[3], 2), "\nSE: ",
                          round(mainres_hit$se[3], 2),
                          '\np = 0.09'),
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Hit Rate",
       title = "A. Hit Rate") + 
  theme_tufte()

resny2 = 
  df_stops_nyc_ts %>% 
  filter(date >= as.Date("2020-04-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-06-15"),
             linetype = 2) + 
  geom_point(aes(x = date, y = arrest),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = arrest, group = por),
              col = 'black',
              se = FALSE,
              size = .6) + 
  annotate('text',
           x = as.Date("2020-07-07"),
           y = .775,
           label = paste0('RDiT Est: ',
                          round(mainres_arr$coef[3], 2), "\nSE: ",
                          round(mainres_arr$se[3], 2),
                            '\np = 0.25'),
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Arrest Rate",
       title = "B. Arrest Rate") + 
  theme_tufte()

resny3 = 
  df_stops_nyc_ts %>% 
  filter(date >= as.Date("2020-04-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-06-15"),
             linetype = 2) + 
  geom_point(aes(x = date, y = rr_bw),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = rr_bw, group = por),
              col = 'black',
              se = FALSE,
              size = .6) + 
  annotate('text',
           x = as.Date("2020-07-07"),
           y = 12,
           label = paste0('RDiT Est: ',
                          round(mainres_rrbw$coef[3], 2), "\nSE: ",
                          round(mainres_rrbw$se[3], 2),
                          '\np < 0.001'),
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Rate Ratio (Black/White)",
       title = "C. Rate Ratio (Black/White)") + 
  theme_tufte()

resny4 = 
  df_stops_nyc_ts %>% 
  filter(date >= as.Date("2020-04-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-06-15"),
             linetype = 2) + 
  geom_point(aes(x = date, y = rr_lw),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = rr_lw, group = por),
              col = 'black',
              se = FALSE,
              size = .6) + 
  annotate('text',
           x = as.Date("2020-07-01"),
           y = 3.5,
           label = paste0('RDiT Est: ',
                          round(mainres_rrlw$coef[3], 2), "\nSE: ",
                          round(mainres_rrlw$se[3], 2),
                          '\np < 0.001'),
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Rate Ratio (Latino/White)",
       title = "D. Rate Ratio (Latino/White)") + 
  theme_tufte()

resny_grob = 
  arrangeGrob(resny1, resny2, resny3, resny4, ncol = 2)
ggsave(plot = resny_grob, filename = 'pics/resny.jpeg', dpi = 600,
       width = 8, height = 4.5)

#### main results: crime ####

load(file = 'data/nyc_crime_cleaned.RData')

nypd_crime = 
  nypd_crime %>% 
  filter(date >= as.Date("2019-01-01")) %>% 
  filter(date <= as.Date("2021-12-31")) %>% 
  group_by(date) %>% 
  summarize(crime_violent = sum(crime_violent, na.rm = TRUE),
            crime_property = sum(crime_property, na.rm = TRUE),
            crime_society = sum(crime_society, na.rm = TRUE)) %>% 
  mutate(por = ifelse(date >= as.Date('2020-06-15'), 1, 0)) %>% 
  mutate(porrv = date - as.Date("2020-06-15")) %>% 
  mutate(blm = ifelse(date >= as.Date('2020-05-28'), 1, 0)) %>% 
  mutate(blmrv = date - as.Date("2020-05-28")) 

# descriptives: por 

porm3 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 1, kernel = 'uni'))
porm2 = with(nypd_crime, rdrobust(x = porrv, y = crime_property, p = 1, kernel = 'uni'))
porm1 = with(nypd_crime, rdrobust(x = porrv, y = crime_society, p = 1, kernel = 'uni'))

pc1 = 
  nypd_crime %>% 
  filter(date >= as.Date("2020-04-15")) %>% 
  filter(date <= as.Date("2020-08-15")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-06-15"),
             linetype = 2) + 
  geom_point(aes(x = date, y = crime_society),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_society, group = por),
              col = 'black',
              se = FALSE,
              size = .6) + 
  labs(x = "Date", y = "Society Crime",
       title = "A. Society Crime (POD)") + 
  annotate('text',
           label = paste0("RDiT Est: ",
                          round(porm1$coef[1], 2), "\nSE: ",
                          round(porm1$se[3], 2),
                          "\np < ", round(porm1$pv[3], 2)),
           x = as.Date("2020-07-15"),
           y = 375,
           family = 'serif',
           size = 2.5) + 
  theme_tufte(base_size = 10)

pc2 = 
  nypd_crime %>% 
  filter(date >= as.Date("2020-04-15")) %>% 
  filter(date <= as.Date("2020-08-15")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-06-15"),
             linetype = 2) + 
  geom_point(aes(x = date, y = crime_property),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_property, group = por),
              col = 'black',
              se = FALSE,
              size = .6) + 
  labs(x = "Date", y = "Property Crime",
       title = "B. Property Crime (POD)") + 
  annotate('text',
           label = paste0("RDiT Est: ",
                          round(porm2$coef[1], 2), "\nSE: ",
                          round(porm2$se[3], 2),
                          "\np = ", round(porm2$pv[3], 2)),
           x = as.Date("2020-07-15"),
           y = 600,
           family = 'serif',
           size = 2.5) + 
  theme_tufte(base_size = 10)

pc3 = nypd_crime %>% 
  filter(date >= as.Date("2020-04-15")) %>% 
  filter(date <= as.Date("2020-08-15")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-06-15"),
             linetype = 2) + 
  geom_point(aes(x = date, y = crime_violent),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_violent, group = por),
              col = 'black',
              se = FALSE,
              size = .6) + 
  labs(x = "Date", y = "Violent Crime",
       title = "C. Violent Crime (POD)") + 
  annotate('text',
           label = paste0("RDiT Est: ",
                          round(porm3$coef[1], 2), "\nSE: ",
                          round(porm3$se[3], 2),
                          "\np = ", round(porm3$pv[3], 2)),
           x = as.Date("2020-07-15"),
           y = 350,
           family = 'serif',
           size = 2.5) + 
  theme_tufte(base_size = 10)

# descriptives: blm protest

blmm3 = 
  with(nypd_crime, rdrobust(x = blmrv, y = crime_violent, p = 1, kernel = 'uni'))
blmm2 = 
  with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 1, kernel = 'uni'))
blmm1 = 
  with(nypd_crime, rdrobust(x = blmrv, y = crime_society, p = 1, kernel = 'uni'))

pc4 = nypd_crime %>% 
  filter(date >= as.Date("2020-03-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2) + 
  geom_point(aes(x = date, y = crime_society),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_society, group = blm),
              col = 'black',
              se = FALSE,
              size = .6) + 
  annotate('text',
           label = paste0("RDiT Est: ",
                          round(blmm1$coef[1], 2), "\nSE: ",
                          round(blmm1$se[3], 2),
                          "\np = ", round(blmm1$pv[3], 2)),
           x = as.Date("2020-07-01"),
           y = 300,
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Society Crime",
       title = "D. Society Crime (BLM)") + 
  theme_tufte(base_size = 10)

pc5 = nypd_crime %>% 
  filter(date >= as.Date("2020-03-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2) + 
  geom_point(aes(x = date, y = crime_property),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_property, group = blm),
              col = 'black',
              se = FALSE,
              size = .6) + 
  annotate('text',
           label = paste0("RDiT Est: ",
                          round(blmm2$coef[1], 2), "\nSE: ",
                          round(blmm2$se[3], 2),
                          "\np = ", round(blmm2$pv[3], 2)),
           x = as.Date("2020-07-01"),
           y = 575,
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Property Crime",
       title = "E. Property Crime (BLM)") + 
  theme_tufte(base_size = 10)

pc6 = nypd_crime %>% 
  filter(date >= as.Date("2020-03-28")) %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2) + 
  geom_point(aes(x = date, y = crime_violent),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = crime_violent, group = blm),
              col = 'black',
              se = F,
              size = .6) + 
  annotate('text',
           label = paste0("RDiT Est: ",
                          round(blmm3$coef[1], 2), "\nSE: ",
                          round(blmm3$se[3], 2),
                          "\np = ", round(blmm3$pv[3], 2)),
           x = as.Date("2020-07-01"),
           y = 350,
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Violent Crime",
       title = "F. Violent Crime (BLM)") + 
  theme_tufte(base_size = 10)

pc_grob = arrangeGrob(pc1, pc2, pc3, pc4, pc5, pc6, ncol = 3)

ggsave(plot = pc_grob, filename = "pics/pcnyc.jpeg", width = 8, height = 4.5)

##### assessing robustness of link between POD and violent crime ####

rdrv1 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 1, kernel = 'uni'))
rdrv2 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 2, kernel = 'uni'))
rdrv3 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 3, kernel = 'uni'))

rdrv4 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 1, kernel = 'tri'))
rdrv5 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 2, kernel = 'tri'))
rdrv6 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 3, kernel = 'tri'))

rdrv7 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 1, kernel = 'epa'))
rdrv8 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 2, kernel = 'epa'))
rdrv9 = with(nypd_crime, rdrobust(x = porrv, y = crime_violent, p = 3, kernel = 'epa'))

crimerobp1 = data.frame(
  
  est = c(rdrv1$coef[1],
    rdrv2$coef[1],
    rdrv3$coef[1],
    rdrv4$coef[1],
    rdrv5$coef[1],
    rdrv6$coef[1],
    rdrv7$coef[1],
    rdrv8$coef[1],
    rdrv9$coef[1]),
  
  se = c(rdrv1$se[3],
    rdrv2$se[3],
    rdrv3$se[3],
    rdrv4$se[3],
    rdrv5$se[3],
    rdrv6$se[3],
    rdrv7$se[3],
    rdrv8$se[3],
    rdrv9$se[3]),
  
  pv = c(rdrv1$pv[3],
    rdrv2$pv[3],
    rdrv3$pv[3],
    rdrv4$pv[3],
    rdrv5$pv[3],
    rdrv6$pv[3],
    rdrv7$pv[3],
    rdrv8$pv[3],
    rdrv9$pv[3]),
  
  spec = factor(c("Uni., 1",
           "Uni., 2",
           "Uni., 3",
           "Tri., 1",
           "Tri., 2",
           "Tri., 3",
           "Epa., 1",
           "Epa., 2",
           "Epa., 3"),
           levels = c("Uni., 1",
                      "Uni., 2",
                      "Uni., 3",
                      "Tri., 1",
                      "Tri., 2",
                      "Tri., 3",
                      "Epa., 1",
                      "Epa., 2",
                      "Epa., 3"))
  
) %>% 
  ggplot() + 
  geom_point(aes(x = spec, y = est),
             size = 2.5) + 
  geom_hline(yintercept = 0) + 
  geom_errorbar(aes(x = spec,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    width = 0)) + 
  labs(x = "Kernel, Polynomial Specification",
       y = "RDiT Coefficient (POD)",
       title = 'A. Effect of POD on Violent Crime') + 
  theme_tufte(base_size = 10)

##### assessing robustness of link between BLM protests and property crime ####

rdrv1 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 1, kernel = 'uni'))
rdrv2 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 2, kernel = 'uni'))
rdrv3 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 3, kernel = 'uni'))

rdrv4 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 1, kernel = 'tri'))
rdrv5 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 2, kernel = 'tri'))
rdrv6 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 3, kernel = 'tri'))

rdrv7 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 1, kernel = 'epa'))
rdrv8 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 2, kernel = 'epa'))
rdrv9 = with(nypd_crime, rdrobust(x = blmrv, y = crime_property, p = 3, kernel = 'epa'))

crimerobp2 = data.frame(
  
  est = c(rdrv1$coef[1],
          rdrv2$coef[1],
          rdrv3$coef[1],
          rdrv4$coef[1],
          rdrv5$coef[1],
          rdrv6$coef[1],
          rdrv7$coef[1],
          rdrv8$coef[1],
          rdrv9$coef[1]),
  
  se = c(rdrv1$se[3],
         rdrv2$se[3],
         rdrv3$se[3],
         rdrv4$se[3],
         rdrv5$se[3],
         rdrv6$se[3],
         rdrv7$se[3],
         rdrv8$se[3],
         rdrv9$se[3]),
  
  pv = c(rdrv1$pv[3],
         rdrv2$pv[3],
         rdrv3$pv[3],
         rdrv4$pv[3],
         rdrv5$pv[3],
         rdrv6$pv[3],
         rdrv7$pv[3],
         rdrv8$pv[3],
         rdrv9$pv[3]),
  
  spec = factor(c("Uni., 1",
                  "Uni., 2",
                  "Uni., 3",
                  "Tri., 1",
                  "Tri., 2",
                  "Tri., 3",
                  "Epa., 1",
                  "Epa., 2",
                  "Epa., 3"),
                levels = c("Uni., 1",
                           "Uni., 2",
                           "Uni., 3",
                           "Tri., 1",
                           "Tri., 2",
                           "Tri., 3",
                           "Epa., 1",
                           "Epa., 2",
                           "Epa., 3"))
  
) %>% 
  ggplot() + 
  geom_point(aes(x = spec, y = est),
             size = 2.5) + 
  geom_hline(yintercept = 0) + 
  geom_errorbar(aes(x = spec,
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    width = 0)) + 
  labs(x = "Kernel, Polynomial Specification",
       y = "RDiT Coefficient (BLM Protest)",
       title = 'B. Effect of BLM Protest on Property Crime') + 
  theme_tufte(base_size = 10)

crimerobp_grob = arrangeGrob(crimerobp1, crimerobp2, ncol = 2)

ggsave(plot = crimerobp_grob,
       filename = 'pics/altspec_crime.jpeg',
       dpi = 600,
       width = 8,
       height = 2.5)

#### main results: assessing 911 calls #### 

load(file = 'data/nyc_911_cleaned.RData')


cfs_nyc_ts = 
  cfs_nyc %>% 
  filter(!(TYP_DESC %in% c('STATION INSPECTION BY TRANSIT BUREAU PERSONNEL',
                        "TRAIN RUN/MOBILE ORDER MAINTENANCE SWEEP",
                        "COMMUNITY TIME"))) %>% 
  mutate(officer_calls = ifelse(TYP_DESC == "VISIBILITY PATROL: DIRECTED", 1, 0)) %>% 
  mutate(civilian_calls = 1 - officer_calls) %>% 
  group_by(date) %>% 
  summarize(officer_calls = sum(officer_calls, na.rm = TRUE),
            civilian_calls = sum(civilian_calls, na.rm = TRUE)) %>% 
  mutate(por = ifelse(date >= as.Date('2020-06-15'), 1, 0)) %>% 
  mutate(porrv = date - as.Date("2020-06-15")) %>% 
  mutate(blm = ifelse(date >= as.Date('2020-05-28'), 1, 0)) %>% 
  mutate(blmrv = date - as.Date("2020-05-28")) 
  
# plotting the difference in officer calls vs civilian calls 

rdout_offcall = 
  with(cfs_nyc_ts, rdrobust(x = blmrv, y = scale(officer_calls),
                            p = 1, kernel = 'uni'))

2600 / mean(cfs_nyc_ts$officer_calls[cfs_nyc_ts$blm == 0])
rdout_offcall$coef[1] * 100

plot_calls_off = cfs_nyc_ts %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  filter(date >= as.Date("2020-03-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2) + 
  geom_point(aes(x = date, y = officer_calls),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = officer_calls, group = blm),
              col = 'black',
              se = FALSE) + 
  annotate('text',
           label = paste0("RDiT Est: ",
                          round(rdout_offcall$coef[1], 2), "\nSE: ",
                          round(rdout_offcall$se[3], 2),
                          "\np < 0.001"),
           x = as.Date("2020-07-01"),
           y = 4500,
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Officer-Initiated Calls",
       title = "A. Officer Calls") + 
  theme_tufte(base_size = 10)


rdout_civcall = 
  with(cfs_nyc_ts, rdrobust(x = blmrv, y = scale(civilian_calls),
                            p = 1, kernel = 'uni'))

960 / mean(cfs_nyc_ts$civilian_calls[cfs_nyc_ts$blm == 0])
rdout_civcall$coef[1] * 100

plot_calls_civ = cfs_nyc_ts %>% 
  filter(date <= as.Date("2020-07-28")) %>% 
  filter(date >= as.Date("2020-03-28")) %>% 
  ggplot() + 
  geom_vline(xintercept = as.Date("2020-05-28"),
             linetype = 2) + 
  geom_point(aes(x = date, y = civilian_calls),
             alpha = .3) + 
  geom_smooth(aes(x = date, y = civilian_calls, group = blm),
              col = 'black',
              se = FALSE) + 
  annotate('text',
           label = paste0("RDiT Est: ",
                          round(rdout_civcall$coef[1], 2), "\nSE: ",
                          round(rdout_civcall$se[3], 2),
                          "\np < 0.001"),
           x = as.Date("2020-07-01"),
           y = 9700,
           family = 'serif',
           size = 2.5) + 
  labs(x = "Date", y = "Civilian-Initiated Calls",
       title = "B. Civilian Calls") + 
  theme_tufte(base_size = 10)


plot_calls_out = arrangeGrob(plot_calls_off, plot_calls_civ, ncol = 2)
ggsave(plot = plot_calls_out, filename = 'pics/callsnyc.jpeg', dpi = 600, width = 8, height = 2.5)
